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Abstract — Image  morphing,  or  image  interpolation  in  the  time 
domain,  deals  with  the  metamorphosis  of  one  image  into  another. 
In  this  paper,  a  new  class  of  image  morphing  algorithms  is  pro¬ 
posed  based  on  the  theory  of  optimal  mass  transport.  The  L 2  mass 
moving  energy  functional  is  modified  by  adding  an  intensity  penal¬ 
izing  term,  in  order  to  reduce  the  undesired  double  exposure  ef¬ 
fect.  It  is  an  intensity-based  approach  and,  thus,  is  parameter  free. 
The  optimal  warping  function  is  computed  using  an  iterative  gra¬ 
dient  descent  approach.  This  proposed  morphing  method  is  also 
extended  to  doubly  connected  domains  using  a  harmonic  parame¬ 
terization  technique,  along  with  finite-element  methods. 

Index  Terms — Image  interpolation,  image  morphing,  image 
warping,  mass  preserving  mapping,  Monge-Kantorovich  flow, 
optimal  transport. 


I.  Introduction 

IMAGE  morphing,  sometimes  referred  to  as  “image  interpo¬ 
lation  in  the  time  domain,”  deals  with  the  metamorphosis  of 
one  image  to  another  [21].  It  is  a  technique  widely  used  in  tele¬ 
vision  commercials,  music  videos  and  motion  pictures.  Image 
morphing  has  also  been  used  for  facial  recognition  [37].  Given  a 
pair  of  images,  the  goal  of  image  morphing  is  to  find  a  sequence 
of  intermediate  images,  such  that  the  first  image  in  the  sequence 
is  equal  to  the  first  given  image  (starting  image)  and  the  last 
image  is  equal  to  the  second  given  image  (ending  image).  The 
process  begins  with  finding  a  reasonable  warping  function  be¬ 
tween  the  two  images,  and  this  warping  function  is  then  used 
to  interpolate  the  position  of  pixels  through  the  in-between  se¬ 
quence.  Finally,  intensity  or  color  interpolation  (i.e.,  cross  dis¬ 
solving)  is  performed  to  generate  the  intermediate  images. 

There  have  been  many  algorithms  proposed  for  image  mor¬ 
phing.  Some  of  the  most  popular  approaches  are  mesh  warping, 
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field  warping,  and  energy-based  warping.  In  mesh  warping  [35], 
features  are  specified  by  a  nonuniform  control  mesh,  and  the 
warping  function  is  usually  generated  by  a  spline  interpolation. 
This  class  of  mesh  warping  algorithm  usually  shows  good  dis¬ 
tortion  behavior,  but  it  has  a  critical  drawback  in  specifying  fea¬ 
tures  since  the  features  on  the  control  mesh  may  have  an  ar¬ 
bitrary  structure.  It  is  also  time  consuming  to  define  the  fea¬ 
ture  correspondence  via  a  user  interface.  In  field  morphing  [2], 
a  pair  of  corresponding  lines  is  specified.  The  mapping  of  a 
point  in  the  vicinity  of  a  line  can  be  determined  by  its  distance 
from  the  line.  In  the  case  of  multiple  line  pairs,  the  warping  of  a 
given  point  is  calculated  by  a  weighted  sum  of  mappings  of  all 
line  pairs.  This  method  is  easy  to  use  to  specify  corresponding 
features.  However,  sometimes  a  part  of  the  image  may  appear 
in  unrelated  regions  in  the  in-between  images  (often  referred 
to  as  “ghosts”).  Energy  minimization-based  warpings  usually 
guarantee  the  one-to-one  mapping  property,  which  prevents  the 
warped  image  from  folding  back  upon  itself.  For  example,  in 
Lee  et  al.’s  work  [21],  points,  polylines,  and  curves  are  sam¬ 
pled  and  reduced  to  a  collection  of  points.  These  points  are  then 
used  to  generate  the  warping  function  by  minimizing  an  energy 
functional.  A  similar  method  has  been  applied  to  facial  mor¬ 
phing  based  on  Navier  elastic  body  spline  [15].  All  of  the  above 
approaches  fall  into  the  landmark-based  category  and  require 
user  inputs  of  the  corresponding  features. 

The  approach  proposed  in  the  present  work  is  intensity  based. 
Here,  one  uses  only  information  given  by  the  images  such  as 
pixel  intensities  to  perform  the  warping.  A  successful  inten¬ 
sity-based  algorithm  can  achieve  automatic  morphing  without 
user  inputs  or  prior  assumptions  on  special  shapes  or  features  of 
the  objects  in  the  image.  See  [31]  for  a  review  of  the  literature 
and  an  extensive  list  of  references  on  this  subject.  We  should 
also  mention  the  nice  work  of  Iwanowski  [17]  in  which  an  image 
morphing  approach  is  proposed  that  combines  morphological 
interpolation  and  linear  filtering  and  does  not  require  control 
points  or  landmarks.  In  this  paper,  we  present  an  automatic  mor¬ 
phing  algorithm  using  a  completely  different  approach  which  is 
based  on  the  theory  of  optimal  mass  transport.  It  is  important  to 
note  that  we  do  not  claim  that  our  method  is  optimal  in  any  sense 
for  morphing,  but  does  throw  some  new  light  on  the  morphing 
problem,  and  seems  reasonable  for  certain  types  of  images  in 
which  there  is  an  elastic  deformation  as  demonstrated  in  some 
of  our  experiments  given  below. 

Optimal  mass  transport  problem  was  first  formulated  by 
Gaspar  Monge  in  1781.  It  concerned  finding  an  optimal  way, 
in  the  sense  of  minimal  transportation  cost,  of  moving  a  pile  of 
soil  from  one  site  to  another.  This  problem  was  given  a  modern 
formulation  by  Kantorovich  in  1948  [19],  and  is  also  known 
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as  the  Monge-Kantorovich  problem  (MKP).  The  optimal  mass 
transport  problem  has  been  extensively  studied  in  the  fields  of 
econometrics,  fluid  dynamics,  automatic  control,  transporta¬ 
tion,  statistical  physics,  shape  optimization,  expert  systems, 
and  meteorology  [26].  Recently,  this  problem  has  been  studied 
within  the  context  of  content-based  image  retrieval  [22],  [27], 
[28].  Pixels  in  an  image  are  divided  into  several  bins  according 
to  their  color  and  spatial  locations.  The  Earth  Mover’s  distance 
(EMD)  is  then  calculated  between  the  bins  of  two  images  and 
used  for  image  retrieval. 

Our  interest  in  MKP  arose  from  our  visualization  work  in 
medical  applications.  For  example,  a  flattened  representation  of 
colon  surface  is  helpful  for  the  detection  of  colon  polys  [13] 
and  a  flattened  representation  of  vessel  surface  is  useful  for 
the  study  of  the  correlation  between  wall  shear  stress  and  the 
development  of  atherosclerosis  [39].  Among  various  flattening 
techniques,  a  mapping  that  preserves  the  area  is  of  special  in¬ 
terests.  In  this  approach,  an  angle  preserving  flattening  is  first 
constructed  through  harmonic  analysis,  and  then  the  optimal 
mass  preserving  mapping  is  applied  in  order  to  correct  area  dis¬ 
tortions.  The  resulting  mapping  is  an  area  preserving  mapping 
with  minimal  distortion.  We  have  also  proposed  an  MKP-based 
image  registration  algorithm  [14],  in  which  a  pseudo-mass  (in¬ 
tensity-weighted  area)  is  preserved.  Our  algorithm  can  handle 
area  preserving  registration  [12]  naturally  by  simply  assigning 
the  pseudo  mass  density  to  unity  on  the  entire  domain  of  the 
image.  Other  approaches  to  using  various  classes  of  diffeomor- 
phisms  for  registration  and  warping  may  be  found  for  example 
in  [25],  [31],  [32],  and  the  references  therein. 

In  this  paper,  we  present  an  improved  approach  for  image 
morphing,  with  special  efforts  taken  to  reduce  the  double  expo¬ 
sure  effect  (also  referred  to  as  the  “fade-in  and  fade-out”  effect). 
An  improved  three-step  gradient  descent  approach  (as  opposed 
to  the  two-step  approach  used  in  [14]  and  [40])  is  employed  here 
for  rectangular  domains.  We  also  explain  how  to  extend  this  ap¬ 
proach  to  irregular  multiconnected  domains. 

We  now  outline  the  contents  of  this  paper.  In  Section  II,  we 
give  a  brief  review  of  the  optimal  mass  transport  problem.  In 
Section  III,  we  present  our  approach  for  image  morphing  be¬ 
tween  two  rectangular  images,  using  the  improved  gradient  de¬ 
scent  method.  Two  types  of  comparison  terms  are  studied  in  this 
section.  In  Section  IV,  we  extend  our  morphing  algorithm  to 
doubly  connected  domains,  based  on  a  harmonic  parameteriza¬ 
tion  technique.  In  Section  V,  we  illustrate  the  proposed  algo¬ 
rithms  using  real  images  as  examples.  Finally,  in  Section  VI, 
we  summarize  the  contributions  of  this  paper  and  discuss  pos¬ 
sible  future  directions.  An  Appendix  is  provided  at  the  end  of 
the  paper  to  give  more  mathematical  details  of  our  methods. 

II.  Monge-Kantorovich  Problem 

We  now  give  a  modern  formulation  of  the  MKP.  Assume  that 
Qo  and  Oi  are  two  domains  in  with  smooth  boundaries, 
each  with  a  positive  density,  ji o  and  (i\,  respectively.  Further, 
we  assume  that  both  domains  contain  same  total  amount  of  mass 

I  =  I  fJL1‘  ^ 

Qq  Oi 


A  mapping  u  :  Qo  — ►  is  called  mass  preserving  (MP),  if  u 

satisfies 

po  =  \Du\pi  o  u.  (2) 

We  denote  this  as  u  e  MP.  Equation  (2)  is  called  the  Jacobian 
equation ,  where  \Du\  denotes  the  determinant  of  the  Jacobian  of 
u,  and  o  denotes  the  composition  of  two  functions.  Equation  (2) 
implies,  for  example,  that  if  a  small  region  in  O0  is  mapped  onto 
a  larger  region  in  Oi ,  there  must  be  a  corresponding  decrease  in 
density  in  order  for  the  mass  to  be  preserved. 

There  exist  infinite  number  of  such  mappings.  A  criterion,  or 
a  metric,  must  be  defined  in  order  to  obtain  an  optimal  map¬ 
ping.  In  this  work,  we  employ  the  Lp  Kantorovich-Wasserstein 
metric  defined  as  follows: 

dp(p o,Pi)p:=  inf  /  \\u(x)  —  x\\p  po(x)dx.  (3) 

u£MP  J 

This  metric  places  a  penalty  on  the  mass  weighted  Lp  distance 
of  each  bit  of  material  moved  by  mapping  u. 

In  this  paper,  we  take  p  =  2.  One  can  show  that  in  this  case 
one  gets  a  unique  minimizer  given  by  the  gradient  of  a  convex 
function;  see  [3],  [11],  [20],  and  the  references  therein.  Note 
that  the  Kantorovich-Wasserstein  metric  defines  the  distance 
between  two  mass  densities,  by  seeking  the  “cheapest”  way  to 
transport  the  mass  from  one  domain  to  another  with  respect  to 
the  metric  (3).  The  optimal  MP  mapping  thus  defined  is  sym¬ 
metric,  the  optimal  mapping  from  CIq  to  Qi  being  the  inverse  of 
the  optimal  mapping  from  Qi  to  Qq. 

III.  Image  Morphing  on  Rectangular  Regions 

In  the  context  of  image  morphing,  we  can  assume  that  the 
mass  density  is  the  image  intensity  and  directly  apply  the 
optimal  mass  transport  algorithm  on  two  related  images  to 
generate  a  deformed  grid.  The  in-between  image  sequence 
can  be  obtained  using  cross  dissolving.  However,  a  mapping 
that  maps  a  small  high  intensity  region  to  a  large  low  intensity 
region  is  not  desirable  since  it  will  cause  double  exposure  effect 
in  the  in-between  images.  The  L 2  Kantorovich-Wasserstein 
metric  imposes  a  penalty  only  on  the  work  spent  on  moving 
mass  from  one  shape  to  another,  but  not  on  the  change  of  mass 
densities  or  image  intensities.  Hence,  we  add  a  comparison 
term  to  the  distance  metric  (3)  to  penalize  the  change  of  inten¬ 
sity  in  the  image.  In  this  section,  we  only  consider  the  image 
morphing  problem  between  two  rectangular  domains,  and  leave 
the  discussion  of  more  general  domains  to  Section  IV. 

The  idea  is  to  minimize  a  functional  of  the  following  form 
over  an  L 2  MP  mappings  u  :  Qq  — >  fij 

Ma[u]  :=  J  C(Io,Iiou)dx-\-a  J \\u(x)  —  x\\2  po(x)dx  (4) 

where  a  £  R  is  a  fixed  positive  number.  The  first  term  controls 
the  “goodness  of  fit”  between  the  intensity  images  Iq  :  Qo  — > 
R  and  I\  o  u  :  Qi  — ►  R,  and  when  /0  and  I\  are  identical, 
this  term  reaches  its  minimum.  The  second  Monge-Kantorovich 
term  controls  the  warping  of  the  map.  The  function  is  the 
mass  density  defined  on  which  could  be  the  same  as  To  or  a 
smoothed  version  of  Jo.  po  could  also  be  any  scalar  field  that  is 
appropriate  for  the  underlying  physical  model.  Similarly,  pi  is 
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the  mass  density  defined  on  Di.  By  adjusting  a ,  we  control  the 
tradeoff  between  minimal  mass  transport  and  minimal  intensity 
change.  It  must  be  pointed  out  that  by  adding  the  comparison 
term  C(7o,  h°  u),  the  energy  functional  (4)  may  have  multiple 
local  minima,  and  the  resulting  optimal  mapping  is  no  longer  a 
curl-free  field  as  in  the  classical  L 2  MKP. 

We  note  that  the  comparison  term  67(70,  I\  o  u)  may  be  taken 
to  be  any  metric  that  measures  the  similarity  between  the  trans¬ 
formed  source  image  and  the  target  image,  e.g.,  the  sum  of 
squared  difference  (SSD),  likelihood  measurement,  correlation 
ratio,  normalized  correlation,  or  mutual  information  (MI)  [34]. 
In  this  paper,  SSD  and  MI  are  adopted  as  the  similarity  mea¬ 
sures,  based  on  the  characteristics  of  our  testing  images. 

If  SSD  is  used  as  the  similarity  measure  [14],  we  are  mini¬ 
mizing  the  following  energy  functional: 

Ma[u\  =  J (Iiou  —  Io)2dx  +  a  J \\u(x)  —  x\\2  /ao(x)dx  (5) 

and  if  MI  is  used  as  the  similarity  measure  [16],  the  energy  func¬ 
tional  has  the  form 

Ma[u\  =  -  j  7, )  loo  Pu  ^JlLdiodh 

.  j.  PIo(lo)Pu  U{ll) 

io  Xi i 

-\-a  /  \\u(x)  —  x\\2  po(x)dx.  (6) 

n0 


There  have  been  a  number  of  algorithms  considered  for  the 
computation  of  an  optimal  transport  map.  For  example,  methods 
have  been  proposed  based  on  Lagrangian  mechanics  that  are 
closely  related  to  the  ideas  in  fluid  dynamics  [3],  and  geometric 
methods  have  also  been  employed  as  in  Culler-Purser  [7] .  Per¬ 
haps  the  most  common  approach  is  to  reduce  the  L 2  optimal 
transport  problem  to  a  linear  programming  problem  [26] .  How¬ 
ever,  a  fundamental  difficulty  of  linear  programming  is  the  com¬ 
putational  complexity  [18]. 

We  propose  to  use  a  different  approach  to  solve  the  modi¬ 
fied  optimal  mass  transport  problem  (4).  The  basic  idea  is  to 
start  from  an  initial  MP  mapping  u°  from  Do  to  Di,  and  based 
on  the  fact  that  the  composition  of  two  MP  mappings  is  still 
an  MP  mapping,  the  second  step  is  to  update  the  mapping  it¬ 
eratively  in  a  gradient  descent  way  using  MP  mapping  s-1  that 
maps  Do  onto  itself  and  satisfies  the  MP  property  to  decrease  the 
pure  L 2  MKP  energy  functional  without  the  comparison  term. 
It  can  be  proved  that  the  curl  of  the  mapping  is  also  decreasing 
in  this  process  [1].  The  resulting  MP  mapping  ft  is  a  curl-free 
field.  Finally,  we  start  from  this  curl-free  field  u  and  use  a  sim¬ 
ilar  gradient  descent  approach  to  minimize  the  modified  energy 
functional  (4).  As  we  mentioned  earlier,  the  functional  (4)  may 
have  multiple  local  minima.  Starting  with  a  curl-free  MP  map¬ 
ping  u  will  make  the  final  MP  mapping  u°°  look  more  natural 
by  putting  most  curl  into  the  dark  region  of  the  images,  as  shown 
in  the  results  described  in  Section  V. 


Note  that  there  is  a  minus  sign  before  the  MI  term,  since  the 
more  similar  the  two  images  are,  the  larger  the  MI  measure  is. 
It  should  also  be  pointed  out  that  the  first  integral  is  taken  over 
the  domain  of  io  x  i  \ ,  where  io  and  ||  are  intensity  levels  in 
the  histograms  of  the  two  images.  The  probability  density  func¬ 
tions  pIo(io )  and  Pu°u(ii)  can  be  estimated  by  the  1-D  non- 
parametric  Parzen-Rozenblatt  density  model  [9] 

pT°  ( i0 )  =  2  f  ip1D  (I0(x)  -  i0)  dx  (7) 

o0 

and 

Pu°u(i l)  “  j'.  /  ^ id(Ii  ( u(x )  -  /, ) dx  (8) 
fto 

where  is  a  1-D  Gaussian  window,  whose  standard  deviation 
can  be  chosen  to  be  10%  of  the  standard  deviation  of  io  and  i\ 
for  (7)  and  (8),  respectively.  V  is  the  area  of  Do,  or  the  number 
of  pixels  inside  Do  f or  the  discrete  case.  From  another  point  of 
view,  this  Parzen- window  method  is  equivalent  to  smoothing  the 
image  histogram  using  a  Gaussian  filter. 

In  a  similar  way,  p 1° ,Jl Ou(io,ii)  can  be  estimated  using  a  2-D 
nonparametric  Parzen-Rozenblatt  density  model  as  follows: 

PuJlOU('io,'ii)  =  y  f  i>2D  (Io(x)  -  io,  h  (u(x))  -  ii)  dx 

o0 

(9) 

where  ^2 d  is  a  2-D  Gaussian  window,  whose  covariance  is  de¬ 
cided  by  the  covariance  matrix  of  the  paired  random  variables 
(io,  i\)-  For  simplicity,  we  write  ^2 d  as  $  in  the  remaining  part 
of  this  paper. 


A.  The  Initial  MP  Mapping  u° 

For  arbitrary  domains,  the  initial  mapping  can  be  solved  using 
the  method  proposed  by  Dacorogna  and  Moser  [8].  Since  we  are 
working  on  more  regular  domains  (specifically,  two  rectangular 
domains  in  2-D  or  two  cubical  domains  in  3-D),  a  simpler  al¬ 
gorithm  is  implemented  here.  To  further  simplify  the  problem, 
we  assume  we  are  working  in  R2  with  Do  =  Di  =  [0,  l]2,  the 
generalization  to  higher  dimensions  being  straightforward.  The 
idea  of  the  proposed  method  is  that  we  solve  a  1-D  mass  trans¬ 
port  problem  in  one  direction  and  then  solve  a  family  of  1-D 
mass  transport  problems  in  the  other  direction.  In  1-D,  the  op¬ 
timal  transport  map  can  then  be  found  by  simple  quadrature. 

We  can  let  the  mass  first  be  transported  along  the  lines  parallel 
to  the  x  axis,  and  then  transported  along  the  lines  parallel  to 
the  y  axis.  Accordingly,  we  assume  that  the  initial  MP  mapping 
has  the  form  u°  =  (a(rr),  b(x ,  y)),  where  function  a  =  a{pc)  is 
defined  by  equation 


a(x)  1 


J  J  fi1(r],p)dr]dp=  J  J  po(rp  p)drjdp.  (10) 


0  0 


0  0 


By  differentiating  (10)  with  respect  to  x,  we  have 
1  1 

a\x)  J  ii\  (a(x),p)  dp  =  J  Ho(x,p)dp.  (11) 


We  may  now  define  function  b  =  b(x,  y)  by  equation 

Kx ,y)  y 

a'  (  ' 

o 


y 

%'(x)  J  Hi(a(x),p)dp  =  J  Ho(x,p)dp.  (12) 


o 
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Since  ay  =  0,  \Du\  =  axby  =  a'(x)by(x,  y),  by  differentiating 
(12)  with  respect  to  y ,  we  get 

a'(x)by(^y)n1  (a(x),b(x,y))  = /j,0(x,y) 

\Du°\/j,i  o  u°  =  no  (13) 

which  is  the  MP  condition.  In  practice,  a(x)  and  b(x.  y)  can 
be  found  using  numerical  integration  techniques.  Given  our  as¬ 
sumption  that  fi o  and  ji\  are  positive  everywhere,  a(x)  is  a 
monotonically  increasing  function.  b(x,y)  is  also  monotoni- 
cally  increasing  with  respect  to  y  from  (12),  given  that  a! (x)  is 
always  positive.  Hence,  there  is  no  space  folding  problem  with 
the  initial  mapping  u°. 

B.  Curl-free  MP  Mapping  u 

Once  an  initial  MP  mapping  u°  is  found,  we  need  to  dis¬ 
cover  the  curl-free  MP  mapping  u  by  finding  the  polar  fac¬ 
torization  [11],  [23]  of  u°.  Rather  than  finding  u  directly,  we 
solve  an  energy  minimizing  problem  iteratively  via  a  gradient 
descent  method.  If  the  energy  functional  is  defined  by  the  L 2 
Kantorovich-Wasserstein  metric,  this  process  is  equivalent  to 
finding  the  polar  factorization  and  is  guaranteed  to  converge  to 
a  global  optimum  [1].  We  give  the  key  steps  of  this  method  in 
this  section,  and  more  mathematical  details  can  be  found  in  [1], 
[14],  and  in  the  Appendix. 

Assume  that  u°  is  updated  by  an  MP  mapping  s  from  Oo  to 
itself,  i.e.,  u  =  u°  o  s_1 .  Based  on  the  facts  that  the  inverse  of  an 
MP  mapping  is  an  MP  mapping  and  the  composition  of  two  MP 
mappings  is  an  MP  mapping,  u  is  also  an  MP  mapping  from  $20 
to  Cl i .  We  take  s  to  be  a  function  of  “time”  (i.e.,  gradient  descent 
parameter),  initially  being  the  identity  map.  In  order  to  maintain 
the  MP  property,  the  update  of  s  should  have  the  following  form: 

Si=(2_c)os  (14) 

V/A)  ) 

for  some  vector  field  (  on  Qq,  with  div(£)  =  0  and  (£,  n)  =  0 
on  dClo,  n  being  the  normal  to  the  boundary  of  Qq.  Accordingly, 
ut  should  satisfy 

uf  = - Du(  (15) 

Mo 

In  the  gradient  descent  approach,  we  take  the  derivative  of  the 
pure  L 2  MKP  energy  functional  with  respect  to  t.  It  can  be 
proved  that  in  2-D,  (  should  satisfy  (  =  V^/,  where  _L  rep¬ 
resents  a  rotation  by  7t/2  counter  clockwise  and  /  is  given  by 

A /  =  —  2  div  [(u  —  id )■*■] 

/=0on  (16) 

Here,  id  indicates  an  identity  mapping.  Considering  (15)  and 
(  =  V^/,  the  updating  equation  for  u  is  given  by 

ut  =  —DuV± A_1div  [(u  —  id )J_]  (17) 

Mo 


C.  inimizer  u°° 

We  use  a  similar  gradient  descent  approach  to  find  a  mini- 
mizer  u°°  to  the  functional  (4).  As  above,  the  idea  is  to  modify 
an  initial  mapping  u°  by  an  MP  mapping  from  TJq  to  itself  in 
order  to  decrease  the  energy  (4). 

More  specifically,  the  updating  equation  is  given  by 

Ut  =  —DuW±A~1div(P±).  (18) 

Mo 


If  SSD  is  used  as  the  comparison  term,  P  has  the  form  of 


Ps  SD  =  f(hou-Io)2Vno  +  —  (Iiou-I0)VIo+2a(u-id) 
Mo  Mo 

(19) 

and  if  MI  is  used  as  the  comparison  term,  P  is  given  by 


P-MI=~y 


1  +  log 


JoJi  ou 


plo  pll  ou 


Viofr) 

Ho(x) 


#a  (Io(x),h  O  U(x)) 

(nIo,I1ou  \ 

1 + ZoAA°v 

=#  I\  o  u(x))  fjpPl 


+  2  a(u  —  id) 


Po(x)  J 


(20) 


where  *  stands  for  2-D  convolution,  and  is  the  derivative  of 
f).  with  respect  to  its  first  variable. 

Once  the  algorithm  converges,  we  have  the  final  warping 
function  u°°.  Then  a  cross-dissolving  method  is  performed 
to  generate  a  sequence  of  in-between  images  7(f),  such  that 
7(0)  =  7o  and  7(1)  =  I\.  It  is  assumed  that  when  time  t  varies 
from  0  to  1,  the  starting  image  7o  continuously  changes  to  the 
ending  image  I\.  We  further  require  that  the  same  transition 
rate  is  applied  to  all  points  on  the  in-between  images  [21]. 
Hence,  the  image  warping  map  at  any  time  t  (t  E  [0, 1])  is 
simply  given  by 


Xf(x)  =  (1  —  t)x  +  tu°°(x)  (21) 

and  the  corresponding  cross-dissolved  image  at  time  t  is  given 
by 


P  {X\x))  =  (1  -  *)/<>(*)  +  th  ( u°°(x )) .  (22) 

7o  and  7i  can  also  be  color  images  and  (22)  can  be  applied 
to  three  color  components  individually.  The  warping  function 
(21)  guarantees  the  continuous  transformation  from  the  source 
image  to  the  target  image,  t  being  the  transition  rate.  One  can 
also  guarantee  that  the  intermediate  frames  are  mass  preserving 
simply  by  shading  the  pixels  in  the  in-between  images  according 

to  i^purVo  o  (x*)-1. 

D.  Implementation 


where  A-1  stands  for  the  solution  to  the  Poisson  equation  (16). 
When  the  algorithm  converges,  we  have  a  curl-free  vector  field 
u ,  which  will  serve  as  the  initial  mapping  for  the  next  step. 


As  a  summary,  our  proposed  algorithm  takes  the  following 
steps. 

1)  Construct  an  initial  MP  mapping  u°  using  (10)  and  (12). 
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2)  Starting  with  u  =  u°,  update  the  MP  mapping  iteratively 
using  (17)  and  obtain  u. 

3)  Starting  with  u  =  u,  update  the  MP  mapping  iteratively 
using  (18). 

4)  Calculate  in-between  images  using  the  cross-dissolving 
method  by  (21)  and  (22). 

An  upwinding  scheme  is  used  when  computing  Du  in  (17), 
and  all  other  spacial  derivatives  (including  V and  div)  are  com¬ 
puted  using  standard  central  differences.  Matlab  solver  poicalc 
is  used  to  solve  the  Poisson  equation  (16).  The  time  step  dt  can 
be  chosen  to  be  less  than 

1  _1 
min  —  (V^A-UivCP-1)) . 

X  ,2  /JjQ 

for  the  nonlocal  flow,  where  i  stands  for  the  component  of  the 
vector.  The  curl-free  mapping  u  is  obtained  as  t  — »  oo.  In 
practice,  the  procedure  iterates  until  the  mean  absolute  curl  is 
sufficiently  small.  More  details  about  the  implementation  of 
Monge-Kantorovich  flows  may  be  found  in  [14]. 

The  numerical  implementation  of  (18)  is  exactly  the  same  as 
that  for  solving  u  using  (17),  except  2 (u  —  id)  is  substituted  by 
P.  The  iteration  stops  when  the  energy  functional  (4)  decreases 
sufficiently  slowly.  The  optimal  mapping  u°°  is  obtained  when 
the  iteration  stops. 

IV.  Image  Morphing  on  Doubly  Connected  Regions 

In  the  previous  section,  we  presented  an  MP  morphing  algo¬ 
rithm  between  two  simply  connected  domains,  or  more  specif¬ 
ically,  two  rectangular  regions.  It  is  assumed  that  mass  is  pre¬ 
served  over  the  entire  domain.  In  some  images,  however,  the 
MP  condition  is  valid  only  between  parts  of  the  two  images 
rather  than  the  entire  domain.  For  example,  if  we  are  given  a 
sequence  of  magnetic  resonance  (MR)  images  of  a  human  heart 
taken  at  different  times  in  the  cardiac  cycle  at  the  same  spatial 
location,  the  MP  condition  is  only  valid  in  the  myocardium  re¬ 
gion,  but  invalid  in  the  heart  ventricles  since  the  volume  of  the 
blood  varies  from  time  to  time.  A  natural  approach  would  be  to 
perform  image  segmentation  first  to  separate  the  myocardium 
and  the  ventricles,  and  then  perform  image  morphing  only  on 
the  myocardial  region  of  interest. 

In  this  section,  we  extend  the  MP  mapping  algorithm  to 
doubly  connected  domains  [41].  The  main  difficulty  comes 
from  the  construction  of  the  initial  MP  mapping.  We  propose 
an  algorithm  that  finds  u°  by  using  harmonic  parameteriza¬ 
tion.  In  this  approach,  the  two  domains  are  first  harmonically 
parameterized,  and  then  u°  can  be  obtained  in  the  harmonic 
coordinate  system.  A  finite-element  method  (FEM)  is  applied 
in  conjunction  with  the  gradient  descent  method  to  account  for 
the  irregularity  of  the  domains. 

A.  Harmonic  Parameterization 

We  first  show  how  to  construct  an  analytic  function  fh  = 
fi  +  if 2  f°r  ^e  harmonic  parameterization  of  a  doubly  con¬ 
nected  domain  (here  i  denotes  a  square  root  of  —1).  Similar 
techniques  have  been  applied  for  colon  surface  visualization 


[13],  tissue  thickness  measurement  [36],  and  defining  orthog¬ 
onal  curves  for  template  matching  [30]. 

Assume  we  have  a  doubly  connected  domain  E,  which  has 
two  boundaries:  the  inner  boundary  a0  and  the  outer  boundary 
a i ,  as  shown  in  Fig.  1. 

The  real  component  of  fh  is  given  by 


A/f  =0 

with  fii&o)  =  0  and  /N^i)  =  1.  (23) 


A  cut  C  is  then  defined  from  g\  to  gq  starting  from  an  arbitrary 
point  x\  on  a i  along  the  negative  gradient  direction  of  ff  and 
meet  the  inner  boundary  go  at  xo.  Cut  C ,  gq  and  g\  form  a  new 
closed  and  oriented  boundary  B  of  domain  E 


7-)  &0  C  (T 1  C 

B  :  xq  — ►  xq  — ►  x\  — ►  x\  — ►  xo 


with  the  constraint  that  E  is  on  the  left  hand  side  of  B.  The 
boundary  condition  of  the  imaginary  component  ff  on  B  is 
given  by 


/a  (0  =  /  =  /  ftT* 


according  to  the  Cauchy-Riemann  equations.  Co  is  any  given 
point  on  B  that  satisfies  (Co)  —  0. 

Finally,  another  Laplace  equation 


A/2  =  0 


(25) 


is  solved  to  give  the  value  of  inside  domain  E. 

Numerically,  we  are  working  on  a  triangulated  domain  and 
the  Laplace  equations  are  solved  by  a  standard  FEM  technique 
as  we  used  for  colon  visualization  [13].  Once  the  analytic  func¬ 
tion  fh  is  obtained,  a  curvilinear  harmonic  polar  coordinate 
system  can  be  defined  by  taking  ff  as  one  coordinate  and  f £  as 
another,  ff  can  be  thought  of  as  a  curvilinear  “radius”  and  as 
the  “angle”.  By  scaling  ff  and  using  a  constant,  we  can  make 
f‘2  run  from  0  to  27r.  Thus,  the  annulus  domain  E  is  mapped 
onto  a  rectangular  region  in  an  angle  preserving  manner.  Fig.  1 1 
shows  such  a  parameterization  on  a  MR  heart  image  with  the 
ventricle  area  excluded. 
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B.  Optimal  Mapping  Between  Doubly  Connected  Domains 

We  follow  exactly  the  same  steps  as  we  used  for  rectangular 
domains  to  find  the  optimal  mapping.  The  first  step  is  to  con¬ 
struct  an  initial  MP  mapping  u°.  By  applying  harmonic  param¬ 
eterization,  (fi0,  t1  o)  is  mapped  onto  a  rectangle  region  (Qq  ?  lA) 
via  conformal  mapping  fh  =  /f  +  if 2  in  an  angle-preserving 
manner.  If  we  define  tA  as 

Mo  =  \Dfhr\o  (26) 

the  mapping  from  O0  to  Oq  is  a  MP  mapping.  Similarly,  domain 
(fill ,  /ij )  is  mapped  onto  another  rectangular  region 
via  gh  =  g\  +  ig\ ,  where  is 

Mi  =  \Dgh\-\i.  (27) 

Since  and  Of  are  both  rectangular  regions  after  the  harmonic 
parameterization,  it  is  easy  to  define  an  MP  mapping  iqnit  be¬ 
tween  them  using  the  method  presented  in  Section  III-A. 

The  whole  process  can  be  illustrated  by  the  diagram  in  the 
equation  shown  at  the  bottom  of  the  page. 

The  resulting  initial  mapping  u°  between  Qq  and  Oi  is  a  com¬ 
position  of  fh ,  and  (gh)  ,  such  that 

U°  =  ( gh )  1  o  Winit  O  fh.  (28) 


(a)  (b) 

Fig.  2.  Two  given  flame  images,  (a)  Starting  image;  (b)  ending  image. 

arrive  at  its  optimal  position  eventually  no  matter  how  it  is  po¬ 
sitioned  initially. 


Obviously,  u°  is  an  MP  mapping  since  fh ,  U[n[ t  and  ( gh )  1  are 
all  MP  mappings. 

Next,  we  evolve  u°  according  to  (17)  and  then  (18)  to  find 
the  optimal  mapping  u°° .  Due  to  the  irregularity  of  the  domains, 
an  FEM  method  is  used  to  solve  the  Poisson  equation  (16)  on 
a  triangulated  surface.  An  upwind  scheme  is  applied  for  com¬ 
puting  Du.  For  all  other  derivatives,  we  use  a  least  mean  squares 
(LMS)  method  to  numerically  calculate  the  first-order  spatial 
derivatives.  For  example,  assume  that  a  given  point  (.To,  yo )  has 
TV  neighbors  (pci,  ^),  i  —  1 ...  TV,  and  a  function  is  defined 
such  that  $(#i,  2/i)  =  for  i  =  0  ...  TV,  it  is  easy  to  see  that 
the  derivatives  of  should  satisfy 

=  (29) 

V  vJ  W-$o/ 

where  A  is  the  position  difference  matrix  given  by 
f  xi-  t0,  yi-yo  \ 

A  =  I  ...  .  (30) 

\xn  -  xo,  yN  ~  yo  J 

When  applying  harmonic  parameterization,  different  selec¬ 
tion  of  cutting  point  x\  may  lead  to  different  initial  mapping 
u°.  However,  the  gradient  descent  method  allows  every  point 
on  the  boundary  to  move  along  the  boundary.  Each  point  will 


V.  Results 

In  this  section,  we  show  the  results  of  applying  the  proposed 
methods  to  various  types  of  imagery.  We  first  test  the  modified 
MP  morphing  with  SSD  and  MI  as  the  comparison  terms  on 
images  of  fire  flames.  The  second  example  is  an  ocean  wave 
video.  By  selecting  several  key  frames  from  the  video,  we  in¬ 
terpolate  the  in-between  frames  using  our  proposed  method  and 
generate  a  continuous  video.  Our  method  seems  well- suited  to 
these  applications  because  the  images  lack  obvious  landmarks, 
require  nonrigid  mappings  for  good  registration.  Further,  the  in¬ 
tensity-based  mass-preservation  constraint  is  reasonable  as  it  re¬ 
flects  physical  reality.  Registration  and  interpolation  of  smoke, 
clouds,  flames  and  fluids  is  an  active  area  of  research  in  the  com¬ 
puter  graphics  and  image  processing  community.  See  [10],  [24], 
and  [33],  for  example,  for  related  methods. 

This  can  serve  as  an  image  compression  method  by  storing 
the  selected  frames  and  generating  the  omitted  frames  on-the-fly 
when  the  video  is  played.  Our  extension  of  the  MP  morphing 
algorithm  to  doubly  connected  domains  is  tested  over  two  MR 
heart  images  acquired  at  the  systolic  and  diastolic  phases  of 
the  cardiac  cycle.  By  applying  the  proposed  method  to  the 
images  with  heart  ventricle  regions  excluded,  we  are  able  to 
generate  a  natural  heart  beat  through  the  cardiac  cycle  using 
only  two  given  images.  The  results  as  images  are  listed  in  this 
paper,  and  the  results  as  videos  can  be  found  on  our  webpage  at 


(^o,^o) 


fh=rt+ifi 


{nlMhi) 


gh=ai+i92 


(Vl  ,Mi) 


'bf'init 
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(1=0,  0.2,  0.4,  0.6,  0.8  and  1) 

Fig.  3.  Morphing  sequence  of  the  flame  imagery  without  comparison  terms. 


(t=0,  0.2,  0.4,  0.6,  0.8  and  1) 

Fig.  4.  Morphing  sequence  of  the  flame  imagery  using  SSD  as  the  comparison  term. 


http  ://w  w  w  .bme .  gatech .  edu/groups/minerva/publications/pa- 
pers/zhu-extra-index.html. 

We  now  explain  the  results  in  more  detail.  The  first  example 
is  image  morphing  between  two  flame  images.  The  starting 
image  [Fig.  2(a)]  and  the  ending  image  [Fig.  2(b)]  are  from 
a  flame  video  sequence  in  the  Artbeats  Digital  Film  Library 
(http://www.artbeats.com).  The  two  images  are  the  24th  and 
the  29th  frames  in  the  video.  Fig.  3  shows  the  result  of  applying 
pure  L 2  MKP  without  a  comparison  term.  Figs.  4  and  5  are 
the  results  of  using  SSD  and  MI  as  the  comparison  term, 
respectively.  Fig.  6  compares  the  results  generated  by  the 
three  different  methods  at  time  t  =  0.5  (the  starting  frame 
is  at  t  =  0  and  the  ending  frame  is  at  /  =  1).  In  the  pure 
L 2  MKP  result  [Fig.  6(a)],  the  double  exposure  effect  is 
obvious.  By  adding  SSD  and  MI  as  the  comparison  term, 
the  double  exposure  effect  has  been  reduced  effectively  as 
shown  in  Fig.  6(b)  and  (c).  The  results  of  SSD  and  MI  are 
very  similar,  however,  since  both  of  the  image  frames  come 
from  the  same  imaging  modality.  MI  has  greater  flexibility 
in  handling  the  case  where  the  two  image  frames  come  from 
different  imaging  modalities,  such  as  PET  and  MRI  in  medical 
imaging,  but  a  better  function  is  needed  to  be  defined  to 
relate  the  image  pixel  intensity  with  the  mass  in  each  of  the 
modalities,  which  is  beyond  the  scope  of  this  paper. 


As  we  have  mentioned  earlier,  adding  a  comparison  term  in 
the  energy  functional  leads  to  the  presence  of  curl  in  the  final 
deformed  grids.  Although  SSD  is  able  to  reduce  the  double  ex¬ 
posure  effect  slightly  better  than  MI,  it  causes  broken  effects  at 
some  edges  between  the  bright  and  the  dark  regions  due  to  the 
high  value  of  curl  remaining  in  these  regions.  Hence,  there  is 
a  tradeoff  between  reducing  the  double  exposure  effect  and  re¬ 
ducing  curl.  In  our  previous  two-step  approach  [14],  [40],  we 
evolved  u  directly  from  u°  according  to  energy  functional  (4). 
As  a  result,  the  remaining  curl  is  mostly  in  the  “bright”  region  of 
the  image  and  it  causes  unnatural  effects  as  a  result.  In  this  work, 
by  starting  with  a  curl-free  mapping  u9  the  algorithm  makes 
corrections  mainly  close  to  the  edges  of  high  intensity  regions 
and  low  intensity  regions.  Since  less  energy  is  required  to  move 
pixels  within  low  density  regions,  the  remaining  curl  is  mostly 
in  the  “dark”  region  of  the  image  and  it  is  almost  unnoticeable. 
These  effects  are  more  obvious  when  observed  dynamically  in 
videos.  Fig.  7  shows  the  deformed  grids  for  these  three  cases. 

In  the  second  example,  a  few  frames  are  selected  from  an 
ocean  wave  video,  with  about  300-ms  intervals  between  two 
consecutive  frames.  These  frames  are  used  as  key  frames  for 
image  morphing.  Our  algorithm  is  applied  to  two  consecutive 
key  frames  to  interpolate  the  motion  of  the  wave.  In  this  ex¬ 
ample,  SSD  is  used  as  the  comparison  term.  The  images  in  the 
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(t=0,  0.2,  0.4,  0.6,  0.8  and  1) 

Fig.  5.  Morphing  sequence  of  th e  flame  imagery  using  MI  as  the  comparison  term. 


(a)  (b)  (c) 

Fig.  6.  Comparing  results  generated  by  pure  L2  MKP,  SSD,  and  MI  at  t  =  0.5.  (a)  Pure  L 2  MKP;  (b)  SSD;  (c)  MI. 


leftmost  column  of  Fig.  8  are  the  key  frames  selected  from  the 
original  video  sequence,  and  an  image  other  then  those  in  the 
leftmost  column  comes  from  the  morphing  result  between  the 
first  image  in  its  row  and  the  first  image  in  the  next  row.  The 
last  key  frame  is  omitted  and  not  shown  here. 

The  third  example  illustrates  image  morphing  over  two 
doubly  connected  domains.  Two  256  x  256  MR  images  of  a 
human’s  heart  were  acquired  using  a  GE  MRI  scanner  at  the 
diastolic  [Fig.  9(a)]  and  the  systolic  [Fig.  9(b)]  phases  of  the 
cardiac  cycle.  Dark  regions  in  Fig.  10(a)  and  (b)  are  two  doubly 
connected  domains,  corresponding  to  the  myocardium  and 
other  tissues  other  than  the  heart  ventricle.  We  also  use  image 
intensity  as  the  mass  density  in  this  case,  given  the  fact  that  MR 
image  intensity  is  a  function  of  photon  density  and,  thus,  related 
to  tissue  mass  density.  The  dark  regions  in  Fig.  10  are  chosen  as 
natural  candidates  to  apply  MP  morphing  to  (in  contrast  to  the 
heart  ventricle  region  in  which  the  change  in  mass  is  too  drastic 
to  satisfy  the  MP  condition).  Harmonic  parameterization  is 
first  performed  over  each  of  the  two  doubly  connected  domains 


(Fig.  11  shows  the  parameterization  over  the  distolic  phase 
image),  and  an  FEM-based  L 2  MKP  is  solved  between  the 
two  domains  to  find  the  correspondence.  In  this  example,  pure 
L 2  MKP  without  comparison  terms  is  adopted.  Fig.  12  shows 
the  final  deformed  grids.  The  in-between  images  can  then  be 
generated  by  cross  dissolving  over  the  entire  domain,  where 
the  deformation  inside  the  heart  ventricle  is  obtained  by  2-D 
spatial  interpolation.  Fig.  13  shows  selected  frames  from  the 
morphing  result. 

VI.  Conclusions  and  Discussion 

In  this  paper,  we  applied  a  modified  Monge-Kantorovich 
flow  to  the  problem  of  image  morphing  by  adding  a  com¬ 
parison  term  to  the  optimal  mass  transport  energy  functional. 
Other  than  the  morphing  examples  presented  in  this  paper, 
this  technique  can  also  be  used  for  the  problem  of  medical 
image  registration,  if  the  underlying  physics  model  satisfies  the 
mass  preserving  condition.  We  are  currently  using  the  image 
intensity  or  a  smoothed  version  of  image  intensity  as  the  mass 
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Fig.  7.  Comparing  deformed  grids  generated  by  pure  L 2  MKP,  SSD,  and  MI.  (a)  Pure  L 2  MKP;  (b)  SSD;  (c)  ML 


density,  which  is  the  simplest  mapping  function  from  intensity 
to  mass.  This  requires  the  initial  and  end  image  to  be  similar 
in  the  sense  that  they  have  the  same  amount  of  mass  in  order 
to  satisfy  the  mass-preserving  condition.  For  images  with  only 
certain  regions  that  satisfy  this  condition,  we  perform  image 
segmentation  first  and  apply  the  morphing  method  only  on  the 
mass-preserving  regions.  This  can  be  done  through  solving  the 
L 2  MKP  between  two  doubly  connected  domains.  If  the  inner 
boundary  is  small  enough,  it  can  be  regarded  as  a  landmark, 
and,  hence,  we  have  solved  the  problem  of  MP  morphing  with 
a  single  landmark.  This  technique  can  be  extended  to  mul- 
ticonnected  domains  (corresponding  to  multiple  landmarks), 
simply  by  dividing  each  domain  into  several  doubly  connected 
domains  to  build  the  initial  mapping  [38]. 

The  L2  mass  moving  penalty  studied  in  this  paper  may  be 
sometimes  too  severe.  It  tends  to  favor  changes  in  density  over 
moving  the  mass  around.  In  fact,  this  is  the  main  reason  that 
causes  the  double  exposure  effects  commonly  seen  in  image 
morphing  applications.  L 1  type  mass  moving  penalty  is  much 
less  severe  for  large  movement,  but  it  may  produce  nonsmooth 
results.  To  solve  this  problem,  we  are  considering  the  use  of 
L1+e  penalty  in  our  future  research,  with  e  being  a  value  be¬ 
tween  0  and  1 . 

The  MP  constraint  can  also  be  combined  with  the  concept  of 
harmonic  mapping  to  provide  a  new  approach  of  MP  diffeomor- 
phisms,  i.e.  to  consider  the  minimization  of  the  Dirichlet  integral 
over  all  MP  maps 

min  [  ||Z>u||2.  (31) 

ueMP  J 

The  underlying  physical  assumption  of  MP  mapping  is 
exactly  the  same  as  that  of  extended  optical  flow  constraint 


(EOFC)  [4].  Developing  new  optical  flow  algorithms  within  the 
framework  of  optimal  mass  preserving  mapping  can  also  be  an 
interesting  future  direction. 

Appendix 

In  this  Appendix,  we  give  mathematical  details  omitted  in  the 
main  sections  of  this  paper.  We  will  mainly  focus  on  the  proof 
of  MP  mapping  properties  and  the  derivation  of  our  improved 
gradient  descent  method. 

A.  Properties  of  Mass  Preserving  (MP)  Mappings 

We  have  stated  in  Section  III-B  that  in  order  to  preserve  the 
MP  condition,  the  evolution  of  s  and  u  must  satisfy 


St  =  (  —  (]  °S 

\IH>  J 

(32) 

ut  = - Dvf 

Mo 

(33) 

for  some  vector  field  (  on  Qq,  with  div(£)  =  0  and  (£,  n)  =  0 
on  <9Q0>  n  being  the  normal  to  the  boundary  of  Qq.  Now  we 
give  the  proof  to  this  statement.  Since  s  is  an  MP  mapping  from 
(Jio?  Mo)  to  itself,  according  to  the  Jacobian  constraint  we  have 
po  =  \Ds\po  o  s.  By  differentiating  it  with  respect  to  time  t ,  we 
get 

0  =  \Ds\tiM)  O  S  +  |Ds|(^0  O  s)t 

0  =  |Ds|  (div(st  o  s_1)  o  s)  no  °  s  +  |Ds|  {(V^o)  °  s,  St) 

0  =  (ii0div(st  o  s-1))  o  s  +  ((V/x o)  o  s,  st) 

0  =  /x0div(st  o  s_1)  +  (V/Uo,  st  o  s-1) 

0  =  div(^os*  °  s_1)-  (34) 


1490 


IEEE  TRANSACTIONS  ON  IMAGE  PROCESSING,  VOL.  16,  NO.  6,  JUNE  2007 


Fig.  8.  Morphing  sequence  of  a  wave  imagery.  Leftmost  column:  selected  key  frames.  Others:  Interpolated  images.  SSD  used  as  the  comparison  term. 


Hence,  st  must  have  the  following  form: 

St=(— c]os  (35) 

V^o  ) 

for  some  vector  field  (  on  with  div(£)  =  0  and  (£,  n)  =  0 
on  dQo,  ft  being  the  normal  to  the  boundary  of  Qq.  Since  u  and 
s  satisfy 

u  o  s  =  u°  (36) 

by  taking  the  derivative  of  (36)  with  respect  to  £,  we  have 

(Du  o  s)st  +  ut  o  s  =  0, 

utos=  —  (Du  o  s)st 
ut  =  —  Dust  O  s_1. 


From  (35),  we  get 

ut  =  -  —Du(.  (37) 

IJ-o 

B.  Gradient  Descent  Methods  for  Pure  L2  MKP 

In  the  following  discussion,  we  intend  to  minimize  a  more 
general  form  of  the  energy  functional 

M[u]  =  /  (u(x)  —  x )  tio(x)dx  (38) 

where  $  is  a  non-negative  C 1  function.  When  <b(;r)  =  \\x\\2, 
it  is  exactly  the  L 2  Kantorovich-Wasser stein  functional.  The 
first  step  is  to  find  an  initial  MP  mapping  u°  as  stated  in  Sec¬ 
tion  III-A.  The  second  step  is  to  minimize  the  energy  functional 


ZHU  et  al. :  IMAGE  MORPHING  TECHNIQUE 


1491 


(a)  (b) 


Fig.  10.  Segmentation  results  of  the  heart  images  in  Fig.  9.  (a)  Diastolic  phase; 
(b)  systolic  phase. 


over  u  =  u°  o  s-1  by  varying  s  over  MP  mappings  from  S20  to 
Qq,  starting  with  s  equal  to  the  identity  map. 

A  change  of  variable  is  applied  here  by  substituting  x  in  (38) 
with  y  —  5_1(rr).  Due  to  the  MP  property  of  s  and  s_1,  the 
following  is  obvious: 

Ho(x)dx  =  Mo  o  s(y)ds(y )  =  Mo  °  s(y)  \Ds(y)\  dy  =  Mo {y)dy 

(39) 

and  also  u(x)  =  uos(y  )  =  u°(y).  Hence,  functional  (38)  equals 
M  =  f  $  ( u°(y )  -  s(y))  Mo (y)dy.  (40) 

By  taking  the  derivative  of  (40)  with  respect  to  t,  we  get 

=  (v®{u0(y)-s(y)),^(y)jiM)(y)dy.  (41) 

Qo 

Then  we  do  another  change  of  variable  by  substituting  y  back 
with  x  =  s(y)  and  get 

=  -  J  (v$  (u(x)  -  x ) ,  Mo (x)st  O  s-\x))  dx.  (42) 

O0 

Clearly,  were  it  not  for  the  constraint  div^o^t  os-1)  =  0,  we 
could  take  yost  o  s-1  =  V$(u(x)  —  x)  to  decrease  energy. 
However,  considering  this  MP  constraint,  yost  o  s_1  should  be 
a  divergence-free  vector  field.  Hence,  we  define 

C  =  Mo(«t  os)-1.  (43) 

V$(n(rr)  —  x)  can  be  decomposed  into  a  curl-free  part  and  a 
divergence-free  part  (Helmholtz  decomposition  [29]),  i.e., 


Fig.  1 1 .  Harmonic  parameterization  of  the  heart  image  in  the  diastolic  phase. 


Fig.  12.  Deformed  grids  over  the  heart  image  in  the  systolic  phase. 


where  div(x)  =  0  and  (%,  n)  =  0  on  dfio-  Then,  (42)  can  be 
rewritten  as 

~Mt=  f{Vw  +  X,0 

O0 

=  J(Vw,C)  +  J  (x,C) 

n0  n0 

=  J  (div(wC)  -  wdiv(C))  +  J (x,  C) 

Qo 

=  J  w((,n)  +  /<*0  =  /<*O  (45) 

where  (  is  chosen  to  be  %,  and  can  be  found  through  the 
Helmholtz  decomposition. 

Now  we  give  the  evolving  equation  for  u  in  the  general  Rd 
case,  as  well  as  in  the  2-D  case  which  has  a  simpler  expression 
and  will  be  used  in  our  algorithm. 

1 )  Gradient  Descent  in  Rrf:  By  taking  the  divergence  of  (44) 
on  both  sides,  it  is  easy  to  see  that  w  is  a  solution  of  the  following 
Neumann-type  boundary  problem 

div  (V4>  (u(x)—x))  =A w 

(Vw,n)  =  (V4>  (u(x)  —  x)  ,n)  on  dQo  (46) 


V4>  (u(x)  —  x)  =  Vip  +  x 


(44) 
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Fig.  13.  Morphing  sequence  of  the  heart  imagery  using  two  given  images  in  Fig.  9. 


and  we  can  set  x  =  S7<fr(u(x)  —  x)  —  Vw.  It  is  then  easily 
seen  that  x  satisfies  the  MP  constraints.  This  PDE  can  be  solved 
using  a  number  of  available  methods,  such  as  the  finite- volume 
method.  Thus,  by  (37),  we  have  the  following  evolution  equation 
for  u: 

ut  — - Du  {V4>  (u(x)—  x)  —  VA-1div  [V4>  {u{x)~  #)]} 

Mo 

(47) 

which  is  a  first-order  nonlocal  scheme  for  ut,  if  we  count  A-1 
as  minus  2  derivatives. 

If  Q(u(x)  —  x )  has  the  form  of  \\u(x)  —  x\\2,  which  is  the  L 2 
MKP,  (47)  has  the  form  of 

2 

ut  = - Du  \u(x)  —  x  —  VA-1div  [u{x)  —  x]}  .  (48) 

Mo 

2)  Gradient  Descent  in  R2:  The  situation  is  simpler  in  the 
R2  case,  due  to  the  fact  that  a  divergence  free  vector  field  ( 
can  be  written  as  (  =  V^h  for  a  scalar  function  M,  where  _L 
represents  rotation  by  7t/2  counter  clockwise,  so  that  S/±h  = 
(—hy,hx).  In  this  case,  (45)  becomes 

-Aft  =  J  <VX/,  Vxft)  =  f  (V/,  Vft)  (49) 

Oo 

where  the  decomposition  of  VT>(,u(rr)  —  x)  is  VT>(,u(rr)  —  x)  = 
ViiJ-f  V1/,  and  we  can  let  h  equal  /.  Hence 

V4>  ('u(rr)  —  —  V/.  (50) 

Considering 

div(V±i6»)  =  div(-ip?/,iprE)  =  -wyx  +  wxy  =  0  (51) 

function  /  can  be  solved  by  the  following  Dirichlet-type 
boundary  problem 

— div  (u(x)  —  rr)^  =  A / 

/  =  0  on  <9^o  (52) 

which  gives  us  the  evolution  equation 

ut  =  —  DuS/±A~1div  fv<I>  (u(x)  —  rr)^  .  (53) 

In  the  L 2  MKP,  (53)  can  be  rewritten  as 

ut  =  — DuV^ A_1div  ((u  —  id)±)  (54) 

Mo 

where  id  is  an  identity  map. 


C.  Gradient  Descent  Methods  for  Energy  Functionals  With  a 
Comparison  Term 

We  now  rewrite  the  energy  functional  as  the  sum  of  two  terms 

M  =  Mi  -f  otM2  (55) 

with  Mi  being  the  comparison  term  penalizing  the  change 
in  intensity.  The  second  term  of  the  energy  functional 
M2  =  J  || u(x)  —  x\\2tiodx  is  exactly  the  L 2  MKP,  which 
has  been  discussed  above.  We  focus  on  the  first  term  Mi,  and 
discuss  the  cases  of  SSD  and  MI,  respectively. 

1)  SSD  as  the  Comparison  Term:  If  SSD  is  employed  as  the 
comparison,  then 

Mi  =  J (Ii  o  u  —  Io)2dx.  (56) 

Before  taking  the  derivative,  we  multiply  Mi  by  po  and  divide 
by  po,  i.e., 

Mi  =  [  —  (Ji  o  u  -  To)2  Mo dx.  (57) 

J  LMo 

By  setting  y  =  5_1(rr)  and  considering  (39):  po {x)dx  = 
Po(y)dy ,  we  get 

Mi  =  J  [h  o  u{x)  -  I0(x)}2  dx 

Mi  =  J  y  [h  °  u(x)  -  Io{x)}2  y,o(x)dx 

Mi  =  [ - [h  o  u°(y )  -  Jo  o  s(y)] 2  y0(y)dy. 

J  Mo  °s(y) 


Then  we  take  the  derivative  of  Mi  with  respect  to  t  and  find  that 


x  Ho (y)dy.  (58) 


By  changing  y  back  to  x,  we  see  that 

-Mit  =  J  [6  °  u(x )  “  -foOO]2  Vflo(x) 

2  d  s  \ 

+~ho(x)  ^  ° ~ Io^  ot  °  dx • 

(59) 
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Adding  the  derivative  of  the  M2  term,  P  is  defined  as 
1  2 

-Pssd  =  —^(Iiou— 7o)2V^oH - (/ion— /o)V/o+2a(w— id) 

Mo  Mo 

(60) 

which  is  (19)  in  Section  III. 

2 )  MI  as  the  Comparison  Term:  If  MI  is  used  as  the  compar¬ 
ison  term 


Mi  =  -  / 


1}.^. 


IqXIi 


P/o(*o)^lOM(*i) 


Now  we  take  the  derivative  of  (61)  respect  to  t  as  in  [16] 


dMt 


d  [  labour-  •  M  Pu  ’/l0“(*o,*i)  J. 

-I 


log 

'nJn 


20  X  u 


-/ 


dp7>’/lOMfa,»i) 

dt 


diodii 


20  X2l 


+ 


/ 


20  X2l 

Considering 


/l  071(2!) 

Pu 


dt 


diodii. 


3Mi 

dt 


/ 


1  +  log 


M 


JoTiou 


(MMi) 


20  X2l 


P/o(*o)p«l0“(*i)J 

,  dPu'hou{iQ,i\) 


Qo 


where  the  same  method  of  multiplying  and  dividing  has  been 
applied.  Then  by  setting  y  =  «s_1(rr)  («s_1  is  an  MP  mapping 
from  O0  onto  itself)  and  noticing  u(x)  =  u°  o  s~*  ,  we  get 


^0,/l°“(«o,*i)  =  2  [  — l—f 
v  J  /.lo  o  s(y ) 

Qo 


(61) 


X  po  (s(l/))  -  *o,  h  (u°(y))  ~  *i))  Po(y)dy.  (66) 
The  derivative  of  p^0,/l0^(«o5  *1)  is  given  as 


<9£ 


-(MMi) 


Qo 


1  fa  (Jo  (s(y))  -  jo,  Ji  (t/°(t/))  -jl) 
k  Mo  0  s(t/) 

#(7)  (s(y)) -*o,  C  (w°(t/))  -«i) 


gf(j/) 

dt 


Mo  °  s(m) 


Po(y)dy 


V/o  (s(j/)) 

Vpo  (s(l/)) 

(67) 


(62) 


where  is  the  partial  derivative  of  f|  with  respect  to  its  first 
component.  Now,  we  do  the  change  of  variable  again  by  substi¬ 
tuting  y  with  s  1(x)  and  derive 

dphJlOU 

Mo-Mi) 


dt 


J  Pu,Il°u(‘io,h)dio  =  Pu°u(ii)  (63) 


it  is  easy  to  see  that  the  last  term  equals  zero;  hence,  (62)  can  be 
written  as  the  following: 


-vf 


fa  (jopr)  -  jp,Ii  (u(x))  -  jj) 
Mo  («) 

f  (Io(x)  -  i0,  h  (u(x))  -  h) 


mx) 


V70(x) 


ds 


—  05  1(x)  \po(x)dx. 


(68) 


Hence,  the  derivative  of  Mi  with  respect  to  t  is  given  by  (69), 
shown  at  the  bottom  of  the  page.  Equation  (69)  is  a  quadruple 
integral,  which  could  be  rewritten  in  a  concise  form  by  using 
convolution 


diodii .  (64)  QM1 


where  pT°,  ph°u,  and  p^:IlOU  can  be  computed  using  (7),  (8), 
and  (9),  respectively.  (d/dt)p{j’IlOU(io,  ii)  can  be  computed  as 
follows.  First,  (9)  can  be  rewritten  as 


dt 


■Kri 

Qo 


1+log 


Jo  ,/i  071 


pl 0  pi  1  ou 


*fa  (Io(x),Il°u(x)) 


»(i)  (  +  Vi>7‘“) 

Vpo(ai) 


x  (/<>(*)  -  «o,  7i  (u(ir)  -  «i))  /j,0(x)dx  (65) 


*f(I0(x),hou(x)) 
^  o  s_1(a;)^  nodx. 


Po(x)  \ 


(70) 


dMi 

dt 


2q  X  2i  Qq 


1  +  log 


P 


Jo,I1OU 


(*o,*i) 


PJ°(*o)p£°^i)J  L 


fa  (Io(x)  -  io,  h  (u(x))  -  ii) 


Mo  09 


V/0(*) 


f  (Iq(x)  ~  ig,h  {u{x))  -  if) 
Mo  0*9 


V^oOO 


d  s 

,  —  o  s_1(a;)  )  /j,0(x)dxdi0di1  (69) 
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Now  it  has  the  form  of  a  double  integral  on  domain,  and  can 
be  combined  with  the  derivative  from  the  M2  term.  Thus,  P  is 


Pmi  = 


1 

V 


1  log 


J0J1OU 


plo  pi  1  ou 


*ipoc(I0(x),I1  o  u(x)) 


xYM$_l1  +  logp 


fojl  ou 


mo  GO 


*ip  {Io{x),h  ou(x)) 


—  2a(u  —  id) 


plo  pl\  ou 

VmoOO 

Mo  60 


(71) 


which  is  (20). 

It  is  easy  to  deduce  the  evolution  equation  for  u,  using  the 
same  approach  for  pure  L 2  MKP.  In  the  Rd  case,  we  find  a 
scalar  function  w  to  be  the  solution  of  the  following  Neumann- 
type  boundary  problem 


div(P)  =  Aw 

(Vw,ft)  =  (P,  n)  on  <9Q0  (72) 


and  we  can  set  x  =  P  —  Vw.  The  evolution  equation  for  u  is 
given  by 

ut  =  — -Dw  {P  —  VA-1div  [V4>  (u(x)  —  #)]}  (73) 

Mo 

with  P  being  either  Pssd  or  Pmi • 

In  the  R2  case,  function  /  can  be  solved  by  the  following 
Dirichlet-type  boundary  problem 


-div(P±)  =A  / 

/  =  0  on  dtt0  (74) 

and  the  evolution  equation  for  u  is  given  by 

ut  =  —DuV±A~1div(P±).  (75) 

Mo 
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